LinearRegressionSlope Function

FUNCTION LinearRegressionSlope &
!
(x, y, nodata, r2) &
!
RESULT (lrs)

IMPLICIT NONE

!Arguments with intent(in):
REAL (KIND = float), INTENT(IN) :: x (:)
REAL (KIND = float), INTENT(IN) :: y (:)
REAL (KIND = float), OPTIONAL, INTENT(IN) :: nodata

!Arguments with intent(out)
REAL (KIND = float), OPTIONAL, INTENT(OUT) :: r2


!Local declarations:
REAL (KIND = float) :: lrs !!slope of linear regression
REAL (KIND = float) :: meanx, meany !!mean of x and y vectors
REAL (KIND = float), ALLOCATABLE :: vx(:), vy(:) !!vectors containing  pairs of valid data
REAL (KIND = float) :: num, den !!numerator and denominator for computing slope
REAL (KIND = float) :: sumx, sumy, sumxy, sumx2, sumy2 !!used to compute r2
INTEGER (KIND = short) :: i, j, count

!---------------------------end of declarations--------------------------------

!check x and y have the same size
IF ( SIZE(x) /= SIZE(y) ) THEN
  CALL Catch ('error', 'Statistics',  &
        'calculate linear regression slope: ', argument = &
        'vectors have different size' ) 
END IF

!remove nodata, keep only pairs of valid data
!count pairs of valid data
IF ( PRESENT (nodata) )  THEN
    count = 0
    DO i = 1, SIZE (x)
      IF ( x (i) /= nodata .AND. &
           y (i) /= nodata ) THEN !found one pair of valid data
           count = count + 1
      END IF
    END DO

    IF (count == 0.) THEN
         lrs = 0.
         CALL Catch ('warning', 'Statistics',  &
            'calculate slope of linear regression: no valid numbers in vectors, ', &
              argument = 'slope set to zero' )
         RETURN
    END IF
END IF

!allocate memory of temporary files and store valid pairs
IF ( PRESENT (nodata) .AND. count > 0 ) THEN
    ALLOCATE ( vx(count), vy(count) )
    j = 1
    DO i = 1, SIZE (x)
      IF ( x (i) /= nodata .AND. &
           y (i) /= nodata ) THEN !found one pair of valid data
           vx(j) = x(i)
           vy(j) = y(i)
           j = j + 1
      END IF
    END DO
END IF

!compute mean of x and y vectors
IF ( PRESENT (nodata) .AND. count > 0 ) THEN
   meanx = GetMeanVector (vx)
   meany = GetMeanVector (vy)
ELSE
   meanx = GetMeanVector (x)
   meany = GetMeanVector (y)
END IF

!cumulate
num = 0.
den = 0.

IF ( PRESENT (nodata) .AND. count > 0 ) THEN
    DO i = 1, count
      num = num + ( vx(i) - meanx ) * ( vy(i) - meany )
      den = den + ( vx(i) - meanx ) ** 2.
    END DO

ELSE
    DO i = 1, SIZE(x)
      num = num + ( x(i) - meanx ) * ( y(i) - meany )
      den = den + ( x(i) - meanx ) ** 2.
    END DO

END IF

!compute slope
IF ( den == 0. ) THEN
  CALL Catch ('warning', 'Statistics',  &
  'calculate slope of linear regression: slope is infinite, ', argument = &
  'slope set to huge number' )
  lrs = HUGE (lrs)
ELSE
  lrs = num / den
END IF

!optional: compute R2 (the square of the correlation coefficient)

IF (PRESENT (r2) ) THEN
    count = SIZE (vx)
    sumx = 0.
    sumy = 0.
    sumxy = 0.
    sumx2 = 0.
    sumy2 = 0.
    DO i = 1, count
        sumx = sumx + vx (i)
        sumy = sumy + vy (i)
        sumxy = sumxy + vx (i) * vy (i)
        sumx2 = sumx2 + ( vx (i) )**2.
        sumy2 = sumy2 + ( vy (i) )**2.
    END DO
   
    r2 = (  count * sumxy - sumx * sumy ) / &
         (  (count * sumx2 - sumx**2. )  * (count * sumy2 - sumy**2. )    )**0.5
    r2 = r2 * r2
    
END IF

IF ( PRESENT (nodata) .AND. count > 0 ) THEN
  DEALLOCATE (vx, vy)
END IF


RETURN
END FUNCTION LinearRegressionSlope